{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "例8.1.1在饲料养鸡增肥的研究中,某研究所提出三种饲料配方:A、是以(鱼粉为主的饲料,4.是以(魄米粉为主的饲料,A,是以直蓿粉为主的饲料.\n",
    "为比较三种饲料的效果,特选24只相似的雏鸡随机均分为三组,每组各喂一种饲料,60天后观察它们的质量.试验结果如表8.1.1所示:\n",
    "                表8.1.1鸡饲料试验数据                          单位:g\n",
    "\n",
    "饲料                         鸡的质量\n",
    "A1     1073     1009    1060    1001    1002    1012    1009    1028\n",
    "A2     1107     1092    990     1109    1090    1074    1122    1001\n",
    "A3     1093     1029    1080    1021    1022    1032    1029    1048\n",
    "\n",
    "本例中,我们要比较的是三种饲料对鸡的增肥作用是否相同.为此,把饲料称为因子,记为A,\n",
    "三种不同的配方称为因子A的三个水平,记为A1,A2,A3,使用配方A下第j只鸡60天后的质量用y。\n",
    "表示,i=1,2,3,j= 1,2,… ,8.我们的目的是比较三种饲料配方下鸡的平均质量是否相等,为此,\n",
    "需要做一些基本假定,把所研究的问题归结为一个统计问题.然后用方差分析的方法进行分析.\n",
    "若相等,可任选一种饲料,特别可选廉价饲料;若不等，应洗增肥效果好的饲料.\n",
    "\"\"\"\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "data": {
      "text/plain": "      1     2     3     4     5     6       7     8\n0  1073  1009  1060  1001  1002  1012  1009.1  1028\n1  1107  1092   990  1109  1090  1074  1122.0  1001\n2  1093  1029  1080  1021  1022  1032  1029.0  1048",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>1</th>\n      <th>2</th>\n      <th>3</th>\n      <th>4</th>\n      <th>5</th>\n      <th>6</th>\n      <th>7</th>\n      <th>8</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>1073</td>\n      <td>1009</td>\n      <td>1060</td>\n      <td>1001</td>\n      <td>1002</td>\n      <td>1012</td>\n      <td>1009.1</td>\n      <td>1028</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>1107</td>\n      <td>1092</td>\n      <td>990</td>\n      <td>1109</td>\n      <td>1090</td>\n      <td>1074</td>\n      <td>1122.0</td>\n      <td>1001</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>1093</td>\n      <td>1029</td>\n      <td>1080</td>\n      <td>1021</td>\n      <td>1022</td>\n      <td>1032</td>\n      <td>1029.0</td>\n      <td>1048</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "productivity_df = pd.read_excel('E:\\Project\\Python\\Statistics\\Third\\chicken.xlsx')\n",
    "productivity_df"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "data": {
      "text/plain": "    factory  time_spent\n0         1      1073.0\n1         1      1107.0\n2         1      1093.0\n3         2      1009.0\n4         2      1092.0\n5         2      1029.0\n6         3      1060.0\n7         3       990.0\n8         3      1080.0\n9         4      1001.0\n10        4      1109.0\n11        4      1021.0\n12        5      1002.0\n13        5      1090.0\n14        5      1022.0\n15        6      1012.0\n16        6      1074.0\n17        6      1032.0\n18        7      1009.1\n19        7      1122.0\n20        7      1029.0\n21        8      1028.0\n22        8      1001.0\n23        8      1048.0",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>factory</th>\n      <th>time_spent</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>1</td>\n      <td>1073.0</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>1</td>\n      <td>1107.0</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>1</td>\n      <td>1093.0</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>2</td>\n      <td>1009.0</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>2</td>\n      <td>1092.0</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>2</td>\n      <td>1029.0</td>\n    </tr>\n    <tr>\n      <th>6</th>\n      <td>3</td>\n      <td>1060.0</td>\n    </tr>\n    <tr>\n      <th>7</th>\n      <td>3</td>\n      <td>990.0</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>3</td>\n      <td>1080.0</td>\n    </tr>\n    <tr>\n      <th>9</th>\n      <td>4</td>\n      <td>1001.0</td>\n    </tr>\n    <tr>\n      <th>10</th>\n      <td>4</td>\n      <td>1109.0</td>\n    </tr>\n    <tr>\n      <th>11</th>\n      <td>4</td>\n      <td>1021.0</td>\n    </tr>\n    <tr>\n      <th>12</th>\n      <td>5</td>\n      <td>1002.0</td>\n    </tr>\n    <tr>\n      <th>13</th>\n      <td>5</td>\n      <td>1090.0</td>\n    </tr>\n    <tr>\n      <th>14</th>\n      <td>5</td>\n      <td>1022.0</td>\n    </tr>\n    <tr>\n      <th>15</th>\n      <td>6</td>\n      <td>1012.0</td>\n    </tr>\n    <tr>\n      <th>16</th>\n      <td>6</td>\n      <td>1074.0</td>\n    </tr>\n    <tr>\n      <th>17</th>\n      <td>6</td>\n      <td>1032.0</td>\n    </tr>\n    <tr>\n      <th>18</th>\n      <td>7</td>\n      <td>1009.1</td>\n    </tr>\n    <tr>\n      <th>19</th>\n      <td>7</td>\n      <td>1122.0</td>\n    </tr>\n    <tr>\n      <th>20</th>\n      <td>7</td>\n      <td>1029.0</td>\n    </tr>\n    <tr>\n      <th>21</th>\n      <td>8</td>\n      <td>1028.0</td>\n    </tr>\n    <tr>\n      <th>22</th>\n      <td>8</td>\n      <td>1001.0</td>\n    </tr>\n    <tr>\n      <th>23</th>\n      <td>8</td>\n      <td>1048.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "productivity_df_long = productivity_df.melt(var_name='factory', value_name='time_spent')\n",
    "productivity_df_long"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "text/plain": "        time_spent                                                           \\\n             count         mean        std     min      25%     50%     75%   \nfactory                                                                       \n1              3.0  1091.000000  17.088007  1073.0  1083.00  1093.0  1100.0   \n2              3.0  1043.333333  43.316663  1009.0  1019.00  1029.0  1060.5   \n3              3.0  1043.333333  47.258156   990.0  1025.00  1060.0  1070.0   \n4              3.0  1043.666667  57.457230  1001.0  1011.00  1021.0  1065.0   \n5              3.0  1038.000000  46.130250  1002.0  1012.00  1022.0  1056.0   \n6              3.0  1039.333333  31.643851  1012.0  1022.00  1032.0  1053.0   \n7              3.0  1053.366667  60.265275  1009.1  1019.05  1029.0  1075.5   \n8              3.0  1025.666667  23.586719  1001.0  1014.50  1028.0  1038.0   \n\n                 \n            max  \nfactory          \n1        1107.0  \n2        1092.0  \n3        1080.0  \n4        1109.0  \n5        1090.0  \n6        1074.0  \n7        1122.0  \n8        1048.0  ",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead tr th {\n        text-align: left;\n    }\n\n    .dataframe thead tr:last-of-type th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr>\n      <th></th>\n      <th colspan=\"8\" halign=\"left\">time_spent</th>\n    </tr>\n    <tr>\n      <th></th>\n      <th>count</th>\n      <th>mean</th>\n      <th>std</th>\n      <th>min</th>\n      <th>25%</th>\n      <th>50%</th>\n      <th>75%</th>\n      <th>max</th>\n    </tr>\n    <tr>\n      <th>factory</th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>1</th>\n      <td>3.0</td>\n      <td>1091.000000</td>\n      <td>17.088007</td>\n      <td>1073.0</td>\n      <td>1083.00</td>\n      <td>1093.0</td>\n      <td>1100.0</td>\n      <td>1107.0</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>3.0</td>\n      <td>1043.333333</td>\n      <td>43.316663</td>\n      <td>1009.0</td>\n      <td>1019.00</td>\n      <td>1029.0</td>\n      <td>1060.5</td>\n      <td>1092.0</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>3.0</td>\n      <td>1043.333333</td>\n      <td>47.258156</td>\n      <td>990.0</td>\n      <td>1025.00</td>\n      <td>1060.0</td>\n      <td>1070.0</td>\n      <td>1080.0</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>3.0</td>\n      <td>1043.666667</td>\n      <td>57.457230</td>\n      <td>1001.0</td>\n      <td>1011.00</td>\n      <td>1021.0</td>\n      <td>1065.0</td>\n      <td>1109.0</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>3.0</td>\n      <td>1038.000000</td>\n      <td>46.130250</td>\n      <td>1002.0</td>\n      <td>1012.00</td>\n      <td>1022.0</td>\n      <td>1056.0</td>\n      <td>1090.0</td>\n    </tr>\n    <tr>\n      <th>6</th>\n      <td>3.0</td>\n      <td>1039.333333</td>\n      <td>31.643851</td>\n      <td>1012.0</td>\n      <td>1022.00</td>\n      <td>1032.0</td>\n      <td>1053.0</td>\n      <td>1074.0</td>\n    </tr>\n    <tr>\n      <th>7</th>\n      <td>3.0</td>\n      <td>1053.366667</td>\n      <td>60.265275</td>\n      <td>1009.1</td>\n      <td>1019.05</td>\n      <td>1029.0</td>\n      <td>1075.5</td>\n      <td>1122.0</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>3.0</td>\n      <td>1025.666667</td>\n      <td>23.586719</td>\n      <td>1001.0</td>\n      <td>1014.50</td>\n      <td>1028.0</td>\n      <td>1038.0</td>\n      <td>1048.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "productivity_df_long.groupby('factory').describe()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "data": {
      "text/plain": "              df       sum_sq      mean_sq        F    PR(>F)\nC(factory)   7.0   7827.18625  1118.169464  0.59554  0.750709\nResidual    16.0  30041.14000  1877.571250      NaN       NaN",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>df</th>\n      <th>sum_sq</th>\n      <th>mean_sq</th>\n      <th>F</th>\n      <th>PR(&gt;F)</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>C(factory)</th>\n      <td>7.0</td>\n      <td>7827.18625</td>\n      <td>1118.169464</td>\n      <td>0.59554</td>\n      <td>0.750709</td>\n    </tr>\n    <tr>\n      <th>Residual</th>\n      <td>16.0</td>\n      <td>30041.14000</td>\n      <td>1877.571250</td>\n      <td>NaN</td>\n      <td>NaN</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "productivity_lm = ols('time_spent~C(factory)', data=productivity_df_long).fit()\n",
    "sm.stats.anova_lm(productivity_lm)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7507087625984523\n"
     ]
    }
   ],
   "source": [
    "data = sm.stats.anova_lm(productivity_lm)\n",
    "p = data.iloc[0,4]\n",
    "print(p)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "因为0.7507087625984523>0.05无法拒绝原假设，\n",
      "即饲料的选择不是重要的影响因素。\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "if p < alpha:\n",
    "    print(\"因为{}<{}拒绝原假设，三种饲料是重要的影响因素。\".format(p,alpha))\n",
    "else:\n",
    "    print(\"因为{}>{}无法拒绝原假设，\\n即饲料的选择不是重要的影响因素。\".format(p,alpha))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}